Moving Finite Size Particles in a Flow: 
A Physical Example for Pitchfork Bifurcations of Tori 
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The motion of small, spherical particles of finite size in fluid flows at low Reynolds numbers is 
described by the strongly nonlinear Maxey-Riley equations. Due to the Stokes drag the particle 
motion is dissipative, giving rise to the possibility of attractors in phase space. We investigate 
the case of an infinite, cellular flow field with time-periodic forcing. The dynamics of this system 
are studied in a part of the parameter space. We focus particularly on the size of the particles 
whose variations are most important in active, physical processes, for example for aggregation and 
fragmentation of particles. Depending on their size the particles will settle on different attractors in 
phase space in the long term limit, corresponding to periodic, quasiperiodic or chaotic motion. One 
of the invariant sets that can be observed in a large part of this parameter region is a quasiperiodic 
motion in form of a torus. We identify some of the bifurcations that these tori undergo, as particle 
size and mass ratio relative to the fiuid are varied. This way we provide a physical example for sub- 
and supercritical pitchfork bifurcations of tori. 

PACS numbers: 05.45.-a, 47.52.-j 
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I. INTRODUCTION 

For a long time there have been intensive studies of the 
chaotic advection of point particles in flows Newer 
studies focus on particles of a finite size, instead of point 
articles without mass. Based on the early work of Basset 
and later Boussinesq [5] and Oseen |(|, equations of 
motion for spherical particles of a finite size in a flow at 
low Reynolds numbers were introduced by Maxey and 
Riley (Maxey-Riley equations). These equations have 
been studied in varying context. In many cases the finite- 
size effects lead to completely different results compared 
to what is observed for ideal tracers In particular 

the results depend on the characteristics of the particles, 
i.e whether they are lighter or heavier than the fluid. 

One of the most important applications of these stud- 
ies are problems where in addition to the dynamics of 
the particles in the flow active processes of the particles 
are included. Such active processes concern the dynam- 
ics of the particles themselves, i.e. they can change their 
physical and chemical properties due to various kinds of 
interaction. Examples are the study of biological pro- 
cesses (growth of plankton populations), chemical pro- 
cesses (chemical reactions in chaotic or turbulent flows), 
and physical processes (aggregation and fragmentation of 
marine aggregates, e. g. s ediment particles in the ocean, 
or cloud formation) |l4l|-(l9j. Using the Maxey-Riley 
equations, the influence of the finite size of particles on 
such active processes as autocatalytic reactions and co- 
alescence of particles has recently been studied in [20| . 
[^l*! . It has been shown that for particular parameters of 
the flow chaotic attractors can occur, leading to filamen- 
tal structures of the particle distribution. 

In this paper the Maxey-Riley equations are studied in 
the case of an infinite, cellular, time dependent flow field. 
Due to the Strokes drag the particle motion is dissipative, 
giving rise to the possibility of attractors in phase space. 



The observed system behavior is highly complex and de- 
pends strongly on the parameter values. Different forms 
of periodic, quasiperiodic and chaotic behavior appear 
in the system, with both single attractors and param- 
eter regions possessing multi-stability. Some particular 
attractors and parameter regions have already been dis- 
cussed by e.g. Maxey Q, Yu et al. [1^ and Nishikawa 
et al. [23] ■ Others have not been studied yet. Because of 
the complexity and the number of parameters that can 
be varied (a total of 6 parameters) only a part of the 
parameter space, namely the variation of two parame- 
ters, is studied in this paper. These two parameters are 
chosen for their relevance in the context of application 
in active, physical processes. In particular we are inter- 
ested in the type of long-term behavior, which occurs for 
different sizes and mass of the particles. Since active pro- 
cesses like aggregation and fragmentation change the size 
of the particles, it is interesting to study the impact of 
variations in the size on the dynamics of the particles. 

For both particles, heavier and lighter than the fluid, 
quasiperiodic solutions in the form of a torus can be ob- 
served in a large parameter region. In this paper we are 
studying some of the changes that these tori undergo, as 
parameters are varied. The main focus is on pitchfork 
bifurcations of tori that appear in this context. 

While there is a large body of work devoted to the 
study of bifurcations of fixed points and periodic orbits 
in various fields of science, there are only a few exam- 
ples from application where bifurcations of tori occur. 
Some general results are known from the mathematical 
literature concerning the bifurcations of tori [23] and par- 
ticularly the break-up of tori and the transition to chaos 
[13], [l^- These theoretical findings are accompanied by 
numerical studies of paradigmatic systems [2^-[2^. Ex- 
amples from physics are mainly focused on the possibility 
of the emergence of tori with three incommensurate fre- 
quencies [l^-jl^ and the break-up of tori with two 
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incommensurate frequencies [33|-[3^- The more basic 
bifurcation of a torus studied here has to our knowledge 
not yet been observed in a physical system and is only 
known from a mathematical point of view (e.g. Sun [st'I). 
For particles moving in a fluid flow two kinds of pitch- 
fork bifurcations of a torus can be obtained in different 
parameter intervals. First we find a subcritical pitch- 
fork bifurcation, where an unstable torus becomes stable 
and gives rise to the emergence of two new unstable tori. 
Second a supercritical pitchfork bifurcation is observed, 
where a stable torus becomes unstable and two new sta- 
ble tori are created. While the first bifurcation happens 
only if the particles are lighter than the fluid, the second 
one occurs for particles heavier than the fluid. In con- 
trast to the usual transition from sub- to supercritical 
which appears on the same bifurcation line, we flnd no 
continuous connection between the two bifurcations. 

The paper is organized as follows. In Section [ll] the 
equations of motions are presented. Based on the Maxey- 
Riley equations the complete equations of motion for 
small spherical particles in a time-dependent flow field 
are recalled and some basic properties of the system 
are described, following the approach of [2l|. Section 
mil gives a general overview of the system behavior de- 
pending on the size of the particles and of the part of 
the parameter space that is of relevance to the pitchfork 
bifurcation presented here. In Section IIV Al the results 
from numerical simulations are described that show a su- 
percritical pitchfork bifurcation of a torus. Section flV Bl 
shows the numerical results for a subcritical pitchfork bi- 
furcation of a torus. Section fVl contains a brief summary. 



II. SYSTEM OF EQUATIONS 

Our investigation is based on the motion of small rigid 
spherical particles in a two dimensional, time periodic 
flow field under the influence of gravity. The flow field 
consists of a regular pattern of vortices, or roll cells, that 
is infinitely extended. This flow field was first introduced 
by Chandrasekhar [ssj as a solution to the Benard prob- 
lem, but has also been used in different context since then 
(e.g. a i). 

The undisturbed flow field is denoted by u{x, t) , where 
X = (xi^X2) is a position in the flow field. The flow is 
incompressible, with constant density p, pressure p and 
dynamical viscosity rj of the fluid. Since we restrict our- 
selves to a two-dimensional flow it can be represented by 
a stream function 

\E'(x, t) — [1 + fc sin(cjt)]^^^ sin(7ra;i/L) • sin(7ra;2/L) . 

TT 

(n.i) 

Uq is the maximum velocity of the flow for a given time 
t, k is the amplitude of the periodic forcing and to is the 
frequency. The size of a vortex is L. 

The velocity field u{x,t) = {ui{x,t),U2{x,t)) of the 
fiuid at the position x is derived from the stream function 



in Eq. ((TLT|) by 

= [Vx^], , 1^1,2 
with * = [0,0,*]. 



(II.2) 



The equations of motion for a particle in the flow de- 
scribed in Eq. (|II.1|1 were derived by Yu et al. A 
brief summary, following the description of Nishikawa et 
al. i^ll], is shown here for completeness. 

The basic equations of motion for the dynamics of a 
small rigid spherical particle of mass nip and radius a in 
an incompressible flow for low Reynolds numbers are the 
Maxey-Riley equations 0] 



Tod 



dV_ 



, Du 
[mp-mFjg + niF— \x(t) 

-lmp-^^{V-u{X{t),t) 



1 



10 

-6TrariY{t) 



1 



d 



'Gna rj I dr- 

/o ^/irv{t ~ t) dT 



f(r)(II.3) 



The position of the particle is denoted by X ~ {Xi,X2) 
and the particle velocity is = {Vi,V2) and Y{t) = 
V{t)-u{X{t),t)-^a^V^u{X{t),t). mp is the mass of the 
displaced fluid, v is the kinematic viscosity and g — gex2 
is the gravitational acceleration, where is the unit 
vector in X2 direction. 

These equations are derived for the conditions 



oWq/v < 1 
{a'/v){Ua/L) « 1 
a/L < 1, 



(11.4) 
(II.5) 
(II.6) 



where Wq is a representative velocity scale for {V — u). 
It is important to distinguish between 



Du du ,^ „, ^ 



(II.7) 



the Lagrangian (substantial) derivative, taken along the 
trajectory of a fluid element and 



du du 
'dt ~ ~dt 



+ {V ■ V)u 



(II., 



the derivative taken along the trajectory of the particle. 
The term 



Du 



is the force from the undisturbed fluid acting on the par- 
ticle at position X{t). 

The buoyancy force is given by 

[nip - mF)g 
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The term 



represents the added mass effect. This term expresses 
the fact that an inertial particle brings a certain amount 
of fluid into motion too. 

The Stokes drag force —6TTariY{t) is proportional to 
the difference between the particle velocity V{t) and the 
fluid velocity u(x, t) . The term 



dr — 

\/T^v{t- 



'dT 



Y{r) 



is the Boussinesq-Basset history term, representing the 
effects of the diffusion of vorticity around the particle. 
This term will be neglected in the following. Manton 
(4^ showed that if fluid inertia effects are included, the 
Basset history term is of less significance. 

The terms with a?'S/^u are the Faxen corrections for 
the nonuniform flow. For the flow field u described by 
Eq. pi .ip , the Faxen corrections are 



(II.9) 



In this case the Faxen corrections only decrease m by a 
small amount, which is proportional to (a/L)^. Because 
of (|II.6|) this term is small and will not affect the quali- 
tative behavior of the system. The Faxen corrections are 
therefore neglected. 

Taking all this into consideration, the equations of mo- 
tion for the inertial particles are reduced to 

mp dV . \^ ^ r-/-^ \ ^''r^ 

+ — j-^ = {mp-mF)g + 6nari[u{X,t)~V\ 

+mp{u ■ V)u + -mp{V ■ V)?! 
3 du 

Introducing dimensionless variables 



and defining 



V 



u 



T^, t* 



^ , (11.10) 
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leads to the dimensionless equations of motion (the as- 
terisks are suppressed for convenience) 

dV(t) - / 1 -\ _^ 3 du 

-^=A[u-V + W]+R(u+-VyVu+-R-. 

(11.12) 



The parameter A represents the effect of the particle iner- 
tia. R is the mass ratio parameter between the fluid and 
the particle. Systems with < | correspond to particles 
heavier than the fluid (aerosols) and i? > | corresponds 
to particles lighter than the fluid (bubbles), i? ^ is 
the aerosol limit, where rUp tends to infinity, compared 
to mi?, i? ^ 2 is the bubble limit, where mp tends to 0, 
compared to m^?. i? = | corresponds to particles which 

are neutrally buoyant. The parameter W is the scaled 
particle settling velocity for still fluid. W is positive for 
bubbles and negative for aerosols. 
Substituting 

LoL 



Uo 



and again suppressing the asterisks yields the dimension- 
less stream function 

\l/(xi(t), a;2(<), i) = — [1 + A:sin(Li;t)] sin(7ra;i) sin(7ra;2) ■ 

(11.13) 

The dimensionless velocity field can then be derived as 



dx2 



u{x{t),t) — I d-fii{t),t) 

\ dxi / 

= [1 + A; sin(cji)] 



(11.14) 



sin(7ra;i) cos(7rx2) 
— cos(7ra;i) sin(7ra;2) 



Combining (|IlT41) with Eq. (ITlT21) results in the full 
equations of motion 



dXi 
dt 

dX2 
dt 



Vi 
V2 



(11.15) 
(11.16) 



dVi 
~dt 



-AVi + A{1 + ksm{Lut)) sin(7rXi) cos(7rX2) 
R 

+— 7r(l + ksm{u;t)){Vi cos(7rXi) cos(7rX2) 
— V2 sin(7rXi) sin(7rX2)) 



i?7r(l + k sm{Lut))'^ sin(7rJY'i) cos(7rXi) 
Ruk cos{u!t) sin(7rXi) cos(7rX2) (11-17) 



dV2 

dt 



-AV2 - A{1 + ksm{ujt)) cos(7rXi) sin(7rX2) 
R 

+— 7r(l + ksin{ujt)){Vi sin(7rXi) sin(7rX2) 
-V2 cos(7rXi) cos(7rX2)) + AW 

+ ^Rn{l + k sm{ujt) f sin(7rX2) cos(7rX2) 
3 

— -i?a;fccos(a;t) cos(7rXi) sin(7rX2) . (11.18) 
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The dynamics of the system Eqs. (jII.15|) - pi.lSI) takes 
place in a five dimensional phase space (two positions, 
two velocities and time). Because of the periodic forcing 
we consider a stroboscopic map M of the system with the 
period T of the forcing. This map projects the dynamics 
onto a four dimensional phase space. 

The dynamics in this four dimensional phase space is 
dissipative. Calculation of the divergence of the flow de- 
fined by Eqs. pi.lSP - (|II.18|) yields that an element 
of the phase space shrinks exponentially with e^^"^*. 
Therefore, attractors can be found in this system. These 
attractors occur in the whole phase space, but in the fol- 
lowing only their projections onto the configuration space 
are shown. 

The dimensionless stream function Eq. has a 

spatial period of 2. This spatial period results in a pe- 
riodic structure of the phase space, with the same pe- 
riod. This periodic structure of the phase space is re- 
fiected in the invariant sets. It can be easily seen that 
for every invariant set of Eqs. (jll.lSp - (jll.lSp there is 
an infinite number of identical invariant sets, repeated 
with the same spatial period 2 throughout the phase 
space. Furthermore, if [Xi{t) , X2{t) jVi{t) ,V2{t)) is an 
invariant set of Eqs. pi.l5[) - (jlLlSp . an identical invari- 
ant set can be found with each of the transformations 
{Xi ^ 2- Xi-Vi ^ -Vi) and {Xi ^ I - Xv^X^ ^ 
1 + X2\ Vi -Vi) and {Xi^l + Xi-X2-^l + X2). It 
is therefore sufficient to consider the map M restricted 
to the unit ceU F = [0, 1] x [0, 1] x M?, with being the 
velocity components. Using the transformations above, 
every invariant set in F can be extended to the whole 
phase space. In the figures shown in this paper, the in- 
variant sets are already extended from the unit cell F to 
the interval F = [0,1] x [0,2] x M^. By extending the 
sets to F it is easier to see what the invariant sets look 
like in the whole phase space, because starting with F 
the sets can be extended in Xi direction by a simple re- 
flection along the axes Xi = 0, 1 and in X2 direction by 
identifying X2 = Q and X2 = 2. Only the projections 
of this interval onto the configuration space {Xi,X2) are 
shown. 

In addition to the parameters R, A and W we intro- 
duce a size class parameter a for the particles. The size 
class parameter a describes how the radius and mass of 
a particle are connected to the radius and mass of a basic 
particle with radius oq and mass toq. We assume that all 
particles are spherical and have a mass and volume that 
are the sum of the masses and volumes of a number a of 
the basic particles. For the radius and the mass of a par- 
ticle of class a this yields = \/a. ■ ao and nia = ct ■ mo 
(with m — nip or m — mp). For the equations of motion 
of the particles this results in a change of the parameter 
A to Aa = a~3A and a change of the parameter W to 
Wa = a^W. 

In the following the parameter values B = 6.4, Q = 
— 1.6, k = 2.72 and uj = n are chosen to allow a compari- 
son with results by e.g. Nishikawa et al. [2l|. The period 
of the forcing is therefore T — 2. The parameters R and 



a are varied. Additionally A and W change accordingly 
when varying R and a. 



III. GENERAL BEHAVIOR 

According to the aim of our study we are interested 
in the qualitative behavior of the nonlinear system de- 
scribed by Eqs. pi.lSp - pi.lSp . Depending on the sys- 
tem parameters we can distinguish different kinds of long- 
term dynamics, such as periodic orbits, quasiperiodic and 
chaotic motion. Our goal is to identify the occurring in- 
variant sets, both stable and unstable ones, and to ana- 
lyze the bifurcations leading to changes in the long-term 
behavior. As already mentioned in the introduction, ag- 
gregation and fragmentation change the size of the par- 
ticles. Even though we do not take these processes ex- 
plicitly into account and deal only with the dynamics of 
passive tracers we focus on the size of the particles as 
the most important bifurcation parameter. We find that 
the dynamics of the particles depends strongly on their 
size, i.e. different invariant sets are obtained for different 
sizes. Additionally we obtain an interesting bifurcation 
that we study in more detail. Before investigating this 
bifurcation we present the overall picture of qualitative 
behavior in some parts of the parameter space. 

Since the system can not be solved analytically, nu- 
merical calculations are required to find the invariant 
sets corresponding to different long-term dynamics. At- 
tractors can be easily found by examining the long-term 
simulations of the particle dynamics. Unstable invariant 
sets are harder to find due to their saddle type character. 
Estimates are made by starting very close to the unstable 
invariant sets and then examining the short-term behav- 
ior. To get a starting point close to the unstable invariant 
set, a refinement procedure has been applied to minimize 
the distance to the unstable invariant set. We follow the 
basic idea of the PIM-triple method which has been de- 
veloped to compute chaotic saddles [1^ . We start with a 
set of initial conditions along a line with fixed X2, V\ and 
V2 but varying X\. Since the unstable invariant set and 
its stable manifolds separate the basins of attraction of 
different stable invariant sets, we integrate all initial con- 
ditions until they reach one of the attractors. Then we 
take the initial condition, say /2, which (i) lies between 
two initial conditions I\ and /s converging to different 
attractors and (ii) possesses the maximum time needed 
to reach the attractor (PIM-proper interior maximum for 
the transient time). This point I2 is then taken as a start- 
ing point for the unstable invariant set. If the transient 
time is not long enough, so that I2 is not close enough to 
the invariant set, then the whole procedure is repeated 
with a line of initial conditions connecting 1\ and /s. 

The qualitative behavior of the system is highly de- 
pendent on the parameter values. In total there are 6 
parameters involved that can be varied. Because of the 
number of different invariant sets that occur for varia- 
tions of each of these parameters no attempt is made 
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FIG. 1: Bifurcation diagram for particles lighter than the 
fluid (mass ratio parameter R = 1). Upper Panel: Projection 
of attractors in one dimension of configuration space. Lower 
Panel: Numerical estimates of the two largest Lyapunov ex- 
ponents Ai (solid line) A2 (dashed line). 



FIG. 2: Bifurcation diagram for particles heavier than the 
fluid (mass ratio parameter R = 0.5). Upper Panel: Pro- 
jection of attractors in one dimension of configuration space. 
Lower Panel: Numerical estimates of the two largest Lya- 
punov exponents Ai (solid line) A2 (dashed line). 



to give a complete description of all the invariant sets 
and bifurcations occurring in the whole parameter space. 
Thus we focus on the variation of a few parameters and 
keep the others fixed. 

As already mentioned the size class parameter a will 
be the main control parameter to be varied. In addition 
we study the behavior for particles with different mass 
ratios compared to the fluid, i.e. for different values of R. 
According to Eq. (jlLlip the variation of R also implies 
a variation of the parameters A and W. 

In Fig. [1] (upper panel) and Fig. [2] (upper panel) bi- 
furcation diagrams (varying a) for two different values 
of R are shown. Fig. [1] is in the bubble region, with 
i? > f, while Fig. m is in the aerosol region, with R < ^. 
Both bifurcation diagrams look rather complicated. For 
different size classes a we observe different long-term dy- 
namics, in particular periodic, quasiperiodic and chaotic 
motion. Varying the size classes we find various bifur- 
cations, indicating transitions between different types of 
motion. Looking only at the projections of the trajecto- 
ries it is not possible to distinguish between quasiperiodic 
and chaotic motion. Such a distinction requires the com- 
putation of the Lyapunov exponents, which are presented 
in Fig. [1] (lower panel) and [5] (lower panel) respectively. 



The computation of the Lyapunov exponents is based on 
the algorithm developed by Shimada and Nagashima [i^l 
which allows for the computation of the whole spectrum 
of Lyapunov exponents separately for each attractor oc- 
curring in the system. It is important to note that the 
figure showing the Lyapunov exponents indicates for each 
parameter value always the largest occurring Lyapunov 
exponents in the system, computed for a set of 30 dif- 
ferent initial conditions. This means that even though 
the Lyapunov exponents are shown as a smooth curve, 
they can correspond to different attractors. Such a visu- 
alization has been chosen to emphasize always the most 
complex motion which can be found for a particular pa- 
rameter value. 

Based on simulations and Lyapunov exponents we can 
immediately identify period doubling, e.g. at i? = 1, a ~ 
4.44 and R = 0.5, a w 6.8, transitions to chaos, e.g at 
R — l,a ^ 4.68 and R — 0.5, a ~ 6.92 and torus bifur- 
cations, e.g. at i? = l,a ~ 0.5 and R = 0.5, a ~ 0.248 
(compare Fig [5] (lower panel) and Fig. [8] (lower panel)). 
Intermittency can be expected at e.g. R = l,a ^ 2.85. 
This is closely related to the emergence of periodic win- 
dows within the chaotic parameter ranges. 

Within the considered parameter range we find regions 
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where only a single attractor occurs and regions with 
multistability where several attractors coexist for the 
given set of parameter values, as i.e. seen for a e [5.2, 5.7] 
for K — 0.5. In the latter regions it depends crucially on 
the initial conditions which of these stable states is re- 
alized. Each of these coexisting attractors has its own 
basin of attraction. For many of the considered param- 
eter ranges these basins have a complexly interwoven 
structure leading to a very high sensitivity of the final 
state to initial conditions (see for example Fig. [3]). In 
the following the basins of attraction shown are com- 
puted as two-dimensional cross-sections of the complete 
four-dimensional basins of attraction for initial condi- 
tions with {Xx,Xi) e [0,1] X [0,2] and (V1,V"2) = (0,0). 




FIG. 3: Black and white dots indicate basins of attraction 
of two different period- 2r orbits for i? = 1, a = 4. The 
period- 2T orbits occur in the whole region [0,2] x [0,2], i.e. 
particles on these orbits move from one unit cell to the next 
(in x\ direction) and then back. The basins of attraction 
display a complexly interwoven structure. To emphasize the 
complicated fractal structure we show in the figure only the 
part [0, 1] X [0, 2] of the configuration space, the part [1,2] x 
[0, 2] is symmetric to the one shown. 

From the physical point of view we note that parti- 
cles of different size behave differently in the same flow, 
i.e. their motion is on distinct attractors. As a conse- 
quence one will find the particles belonging to different 
size classes at different locations in the flow. We also 
note that bubbles and aerosols of the same size exhibit 
different behavior as well. 

From all the occurring bifurcations we choose the torus 
pitchfork bifurcation for a detailed analysis, since to our 
knowledge it has not been discovered in this context of 
application so far. Therefore, we restrict ourselves to 
the parameter region a e [0, 1] and R E [0.4, 1.1]. This 
parameter interval covers parts of the bubble region and 
of the aerosol region in parameter space, where tori and 
their bifurcations can be observed. Both a supercritical 
and a subcritical pitchfork bifurcation of these tori can 
be observed. Note that we study the stroboscopic map 
M in which the quasiperiodic motion on a torus in the 
original flow appears as an invariant curve. 



In the supercritical pitchfork bifurcation a stable torus 
looses its stability and two new stable tori are cre- 
ated. Such a bifurcation has been studied in the quasi- 
periodically forced circle map [i^. There, this bifur- 
cation simply occurs because the quasiperiodic forcing 
turns the normal pitchfork bifurcation of fixed points 
(periodic orbits of period 1) into a pitchfork bifurcation 
of invariant curves in the map, corresponding to tori in 
flows. 
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FIG. 4: Part of the parameter space spanned by size class pa- 
rameter a and the mass ratio parameter R, where the pitch- 
fork bifurcation of a torus occurs. Solid lines indicate a pitch- 
fork bifurcation (supercritical in the aerosol region, subcritical 
in the bubble region). The dashed line marks a torus bifur- 
cation at the border between aerosol and bubble region (at 
i? = 2/3). 

The subcritical pitchfork bifurcation of tori follows a 
different scenario. An unstable torus becomes stable 
and two unstable tori are created which coexist with 
the stable torus. Though proofs of existence for both 
of these bifurcations are known from mathematical lit- 
erature (Broer et al. [1^, Sun [sj) they have not been 
analyzed in detail in the context of a physical application. 

Let us now discuss the occurrence of the pitchfork bi- 
furcation of tori in the two dimensional parameter space 
spanned by R and a. First we note that the pitchfork bi- 
furcation of tori occurs in the bubble region as well as in 
the aerosol region, but there is a fundamental difference 
between them. In the aerosol region (i? < |) the pitch- 
fork bifurcation of a torus is always supercritical, giving 
rise to the emergence of two stable tori. The pitchfork 
bifurcation of a torus in the bubble region (i? > |) is 
always subcritical, with two unstable tori branching off. 

There is a gap between the two bifurcation curves 
marked by the line i? = |, where the particles are neu- 
trally buoyant, so that there is no continuous transition 
from sub- to supercritical as known for many other bi- 
furcations. The two stable tori that are found in the 
aerosol region after the pitchfork bifurcation disappear 
in a Neimark-Sacker bifurcation at 7? = | when R is 
increased (dashed line in Fig. 2]). There each of the sta- 
ble tori merges with an unstable periodic orbit. At that 
point the tori disappear and the periodic orbits become 
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stable. At the same point the spatial symmetry of the 
unstable torus is reversed (compare Fig. [S]and[ni). 

The different cases occurring in this parameter region 
are described in detail in the following section. Both 
the stable and the unstable invariant sets are estimated 
numerically and their projections onto the configuration 
space are shown. 



IV. PITCHFORK BIFURCATION OF TORI 



A. Supercritical pitchfork 
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FIG. 5: Bifurcation diagram and numerical estimates of the 
two largest Lyapunov exponents for R = 0.5, in the area of 
the supercritical pitchfork bifurcation. The arrow indicates 
the bifurcation point at Ucrit = 0.2475. Beyond the bifur- 
cation point some larger periodic windows, corresponding to 
phase-lockings on the tori, can be easily seen. These periodic 
windows appear where the largest Lyapunov exponent drops 
suddenly below 0, e.g at a ~ 0.268. This is also visible in the 
bifurcation diagram (upper panel). 




stable torus 
unstable torus 
unstable fixed 
point 



X 



FIG. 6: Supercritical pitchfork bifurcation of a torus at i? = 
0.5 (solid line - stable torus, dashed line - unstable torus, 
cross - unstable periodic orbit). Upper panel: Before the 
bifurcation (a — 0.2). Lower panel: After the bifurcation 
(a = 0.28). 




FIG. 7: Basins of attraction for the two stable tori for R = 0.5, 
a = 0.325. The basins of attraction are separated by an 
unstable torus. 



In this section the aerosol region of the part of the 
parameter space shown in Fig. 3] is analyzed. In this 
parameter region a supercritical pitchfork bifurcation of 
a torus is observed. The corresponding enlarged bifur- 
cation diagram for R — 0.5 is presented in Fig. [51 A 
solid line in Xi direction for a fixed a indicates a stable 
torus or chaos, while individual points indicate a stable 



periodic invariant set. 

For small values of a the only stable attractor in phase 
space is a single torus. In addition there exist two un- 
stable periodic orbits of period T. An illustration of this 
can be seen in Fig. [5] (upper panel) for a = 0.2. All 
particles end up on the torus. Note that Fig. [S] shows 
the Pointcare-section, therefore the torus appears as an 
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invariant curve and the periodic orbits as fixed points. 
In this aerosol region the torus, and therefore the motion 
of the particles, tends to go through the local downflow 
regions. 

When a is increased, a bifurcation of the stable torus 
takes place at acrit ~ 0.2475. The stable torus splits into 
two stable tori, with an unstable torus separating the 
stable tori (see Fig. [6] (lower panel) for a = 0.28). This 
is a supercritical pitchfork bifurcation of a torus. This 
change can also be seen in the Lyapunov exponents (Fig. 
[5] (lower panel)), where the smaller Lyapunov exponent 
also becomes zero at the bifurcation point. 

Beyond the pitchfork bifurcation it depends on the ini- 
tial condition to which torus the particles converge to. 
The unstable torus can also be seen in the basins of 
attraction of the stable tori, as it makes up the basin- 
boundary together with its stable manifolds (Fig. [7]). 

Beyond the pitchfork bifurcation the behavior gets 
more complicated, as parameter regions with various 
phase-lockings on the unstable and stable tori are en- 
countered. Some of the larger regions with phase-locking 
on the stable tori can be seen in the bifurcation diagram 
Fig. [5l for example at a « 0.268 and 0.305, where the 
zero Lyapunov exponent becomes negative. The param- 
eter regions with phase-locking are separated by large 
parameter regions of stable tori without phase-locking. 

For lower values oi R {R < 0.4), the behavior of the 
system is also much more involved and will not be ex- 
amined in detail here. In this region, a number of bifur- 
cations, in particular phase-lockings, occur on the torus 
before the pitchfork bifurcation. This gives rise to a com- 
plex series of bifurcations, that will be published else- 
where. 

If R is increased and comes closer to the bubble re- 
gion (i? ^ |), another bifurcation of the stable tori can 
be observed. The two stable tori that were created in 
the pitchfork bifurcation move closer to the unstable pe- 
riodic orbits and finally disappear in a Neimark-Sacker 
bifurcation (dashed line in Fig. [4]). There the periodic 
orbits become stable. All particles now move to one of 
these two periodic orbits, with their basins of attraction 
separated by the remaining unstable torus. 



B. Subcritical Pitchfork 

Now we analyze the bifurcation scenarios in the other 
part of the parameter space, namely the bubble region 
with i? > |. The dynamics is dominated by two stable 
fixed points, whose basins of attraction are separated by 
an unstable torus. However, the symmetry of the torus 
is reversed, with the torus now going through the regions 
with local upflow. 

As shown in Fig. 3] we find again a pitchfork bifurca- 
tion line for a torus, but in the bubble region this pitch- 
fork bifurcation is subcritical. 

Let us discuss the bifurcation diagram for R = 1 (Fig. 
[5]) in the vicinity of the bifurcation. For small a values 
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FIG. 8: Bifurcation diagram and Lyapunov exponents for 7? ; 
1, in the area of the subcritical pitchfork bifurcation. 



the long-term dynamics yields two periodic orbits of pe- 
riod T. In addition the dynamics show the existence of 
an unstable torus (see Fig. [5] (upper panel)). The basins 
of attraction of the periodic orbits are separated by the 
unstable torus, resulting in a smooth basin boundary. 

At acrit ~ 0.4996 the unstable torus undergoes a bi- 
furcation. The unstable torus becomes stable, with two 
new unstable tori branching off (see Fig. [5] (lower panel) 
for a — 0.51). This is a subcritical pitchfork bifurca- 
tion of a torus. In the bifurcation diagram Fig. [5] this 
transition is manifested by the abrupt appearance of the 
quasiperiodic motion indicated by the black area beyond 
the critical parameter value. The change can also be seen 
in the Lyapunov exponents, where the largest exponent 
Ai now is equal to zero (compare Fig. [T] (lower panel)). 

There are now three possible states for the long-term 
dynamics, with the actual result for a particle depending 
on the initial condition. Once again, the unstable tori 
form the borders of the basins of attraction. In this case, 
the unstable tori separate the basins of attraction of the 
stable torus from the basins of attraction of the two peri- 
odic orbits (Fig. [T0|) . Here, the basin boundary is again 
smooth. 

For different values of a, the system shows a number 
of other interesting types of behavior, but the analysis 
of these is beyond the scope of this paper and is left for 
further studies. Some examples include a series of phase- 
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FIG. 9: Subcritical pitchfork bifurcation of a torus at i? = 1 
(solid line - stable torus, dashed line - unstable torus, bullet 
- stable periodic orbit). Upper Panel: Before the bifurcation 
(a = 0.4992). Lower panel: After the bifurcation [a = 0.51). 



X^ 
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FIG. 10: Basins of attraction of the stable torus (white) and 
the fixed points (black) for R = \, a = 0.51. The basins of 
attraction are separated by the two unstable tori. 



lockings on the unstable torus and subsequent pitchfork 
bifurcations of the resulting periodic orbits in the inter- 
val a G [0.41,0.45], emergence of other stable periodic 
solutions (a > 0.6) and a break-up of the stable torus at 
a = 1 (see ^ij ) that gives rise to a chaotic attractor with 
fractal basin boundaries. 



V. SUMMARY 

We have investigated the dynamics of passive finite- 
size particles in a fluid flow at low Reynolds numbers. It 
was found that the dynamics of such particles can vary 
greatly with their size and also with their mass relative 
to the fluid. For different parameter regions, periodic, 
quasiperiodic and chaotic motion was observed. When 
considering sets of particles belonging to different size 
classes and advected by the same flow, finite-size particles 
will settle on different attractors depending on their size. 
That means that in systems where aggregation and frag- 
mentation processes take place the particles of different 
size will concentrate in different areas of the configuration 
space. In the case when the time span between the ag- 
gregation and fragmentation events is long compared to 
the time needed to approach the attractor, the particles 
will be distributed according to their respective attrac- 
tors. But in most cases one has to expect that aggrega- 
tion and fragmentation will happen more frequently, so 
that the transient dynamics will blur the separation of 
different size particles in different attractors. 

For small size classes we find a parameter range where 
the dynamics is dominated by quasiperiodic attractors 
for particles lighter and heavier than the fluid. This 
quasiperiodic motion in form of a torus has been stud- 
ied in more detail. The torus undergoes both a sub- and 
a supercritical pitchfork bifurcation, depending on the 
mass of the particles relative to the fluid. For particles 
lighter than the fluid the bifurcation is always subcritical, 
while for particles heavier than the fluid the bifurcation 
is supercritical. No continuous transition between these 
two bifurcations exists in this system, instead the transi- 
tion between the different states happens via a Neimark- 
Sacker bifurcation where the two stable tori emerging in 
the supercritical pitchfork-bifurcation disappear. This 
Neimark-S acker bifurcation forms a boundary at the line 
i? = I of neutrally buoyant particles, separating the dif- 
ferent invariant sets observed for both kinds of pitchfork 
bifurcations. 

At this line of neutrally buoyant particles the symme- 
try of the unstable torus is reversed (compare Fig. [5] 
(upper panel) and Fig. [S] (upper panel)). This change 
in the symmetry can be easily understood from a phys- 
ical point of view. Particles lighter than the fluid get 
suspended at the stable fixed points for long times, but 
for short times (if they start close to the unstable torus) 
move upwards through the flow. The torus therefore lies 
in the regions with local upflow. Particles heavier than 
the fluid sink through the fluid, the torus lies in regions 
with local downflow. At the line i? = | particles neither 
rise nor sink through the flow, instead they collect in the 
region with only horizontal fluid velocity, i.e. here the 
tori become a horizontal line at the border of each unit 
ceU {x2 = 0,1,2). 

For R very close to | the line of supercritical pitchfork 
and the line of Neimark-Sacker bifurcations seem to touch 
tangentially in a bifurcation of higher codimension. But 
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the appearance of higher codimcnsion bifurcations is be- 
yond the scope of this paper. AdditionaUy, for low and 
high values of R {R < 0.4, i? > 1.1) other bifurcations 
occur which are related to a phase-locking on the torus 
under consideration, giving rise to a complex network of 
bifurcation lines in parameter space. 

As a result, a system of finite-size particles moving in 
a fluid flow, provides an example of a sub- and supercrit- 
ical pitchfork bifurcation of a torus in a realistic physical 
system. 
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